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RESUMEN 



El trabajo descrito a continuación, fue realizado con el objetivo de 
encontrar una solución más precisa a los problemas de vigas cuando se utiliza el 
método de los elementos finitos, como también utilizar el programa MathCAD para 
implementar esta solución. Por lo que se propone como primera instancia la 
aproximación de un polinomio de grado cinco para la función de desplazamiento 
transversal de la viga. A lo largo del trabajo, se mostrará el desarrollo detallado de 
la teoría y luego ejemplos de aplicación como demostración de las ventajas para 
su uso. 



ABSTRACT 



The work described below was conducted in order to find more precise 
solutions to the typical beam problems when finite element method is used, as well 
as using the MathCAD program to implement this solution. It is proposed as a first 
approximation for the transversal displacement of the beam, a polynomial with 
degree five for the function of transverse displacement of the beam. Throughout 
the work, detailed development of the theory will be presented, as well as 
application problems, showing the advantage of the use of the given method. 

1 INTRODUCCIÓN 



Los textos que abordan la teoría de elementos finitos en el apartado 
correspondiente a vigas, utiliza una función cúbica para describir el 
desplazamiento. Esta función genera resultados exactos cuando las cargas son 
puntuales, pero cuando la carga es distribuida existe una incompatibilidad en el 
momento ya que al ser este igual a la segunda derivada del desplazamiento se 
generara una representación lineal del mismo. Esta dificultad es solventada 
conforme se aumentan el número de nodos o lo que es lo mismo discretizando la 
viga, conforme se incremente el número de elementos finitos los resultados 
obtenidos serán más cercanos a los resultados que predice la teoría de vigas. 
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La alternativa que se propone es incrementar el orden de la función de 
desplazamiento mediante la inclusión de un nodo interno. Por lo tanto la propuesta 
es desarrollar una función de alto orden de modo que con apenas uno o dos 
elementos finitos dependiendo del problema, se alcance mejor exactitud que con 
el método tradicional de solución (polinomio cúbico). Los pasos a seguir serán 
básicamente proponer una nueva función de desplazamiento incrementando el 
grado del polinomio, derivar las funciones de interpolación N t , desarrollando una 
matriz de rigidez que incrementará la exactitud de cualquier problema de vigas. 
Adicionalmente estos problemas se pueden resolver ventajosamente con las 
herramientas de cálculo simbólico del software MathCAD proporcionando un 
contexto pedagógico muy amplio para la explicación conceptual del método de los 
elementos finitos. 



2 PRINCIPIOS TEÓRICOS 



En el presente numeral se revisara la metodología para obtener la 
formulación del elemento finito del elemento viga en su forma clásica. Una viga es 
un elemento largo y delgado, sujeto a carga transversal que produce significativos 
efectos de flexión, esta flexión es medida como un desplazamiento transversal y 
una rotación, por lo tanto los grados de libertad por nodo son dos. Para un 
elemento finito de longitud L, se tendrán f^m^ái^i, ver figura 1 . 




Figura 1 Deformación de una viga 

Dónde: f ly , corresponden a fuerza y momento local del nodo 1 . 

corresponden a desplazamiento y rotación del nodo 1. 



Se tiene por lo tanto que la ecuación diferencial que gobierna una viga, se 
obtiene del desarrollo de la viga de Euler-Bernoulli, cuyo principio es que las 
secciones planas, inicialmente perpendiculares al eje de la viga, seguirán siendo 
perpendiculares al eje de la viga, aun cuando este se haya curvado, ver Figura 2b. 
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Figura 2 a) viga sin deformar b) Viga deformada debido a carga aplicada c) Elemento 

diferencial de una viga 1 . 



En primer lugar se efectúa el análisis de equilibrio del elemento diferencial, 
mostrado en la figura 2c. 




£ Fy = 0; V - V - dV - w(x) * dx = 0; w(x) = - ^ (1 ) 

I M = 0; - V * dx + dM + w(x) * dx * (y) = 0 o V = ^ (2) 

Donde V es la fuerza cortante 
M es el momento flector 
w(x) es la carga distribuida 

Recurriendo a la teoría de vigas, se determina que la ecuación constitutiva 
de la viga que relaciona solicitaciones con desplazamiento estará definida por: 



Donde k es la curvatura y p es el radio de curvatura 
E es el módulo de elasticidad 
I es el momento de inercia 



Por otro lado la distorsión angular según la figura (3) está dado por: 

dv(x) (4) 
dx 



1 Logan, Daryl; A FIRST COURSE IN THE FINITE ELEMENT METHOD,[UN 
PRIMER CURSO EN EL MÉTODO DE ELEMENTOS FINITOS] 
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Donde: v(x) es el desplazamiento transversal de la viga y 4>, la rotación. 
Adicionalmente se conoce que la curvatura se la define según (5): 

d 2 v(x) 



K = 



dx 2 



(5) 



Si para pequeña curvatura se tiene la ecuación (6): 

_ d 2 v(x) _ M 
K_ dx 2 El 

Por lo tanto el Momento es igual a (7): 

dx 2 

La fuerza cortante sería la derivada del momento según (8): 



V= — = El d3y(x) 

dx dx 3 



Y la carga distribuida w(x) la derivada del cortante, según (9): 

dV d 4 v(x) 
w(x) = - — = -E I — -V 
v ' dx dx 4 



(6) 



(7) 



(8) 



(9) 



Se debe escoger por lo tanto un polinomio con cuatro constantes ya que 
este es el número de grados de libertad del elemento y de grado 3 que satisfaga la 
ecuación diferencial. 



v(x) = al x 3 + a2 x 2 + a3 x + a4 



(10) 



Con las condiciones de frontera apropiadas Fig. 1, se determinan las 
expresiones para las constantes a¡ de (10) 

v(x)= ¿(2x 3 - 3Lx 2 + L 3 ) dly + ¿(x 3 L - 2L 2 x 2 + xL 3 ) c|>l + ^(-2x 3 + 3Lx 2 ) d2y + 

^(x 3 L-L 2 x 2 )c|>2 (11) 



Esta expresión se la reemplaza en (7) y (8), para ambos extremos de la 
viga y se obtiene finalmente la formulación del elemento finito que se va a 
denominar formulación clásica: 
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OBTENCIÓN DE LAS FUNCIONES DE INTERPOLACIÓN N t 



Se propone ahora utilizar un elemento de alto orden con un nodo interno. 
Cuando se incrementa el orden del polinomio se dice que el método de 
refinamiento es "p-refinement". Se sabe por lo tanto que el elemento finito tendrá 6 
grados de libertad, 2 por nodo, sujeto a fuerzas transversales y momentos (Figura 
4) y que el nodo interno no se conecta con ningún elemento. 



9.6 





¿2- rt >¡ 




^f2y,d2y 






í ^ — ^ 



f1y, d1y 



f3y, d3y 



Figura 4 Elemento finito propuesto 



Si se tienen 6 grados de libertad el polinomio de interpolación deberá ser 
de grado 5 de la forma siguiente: 

v(x) = al x 5 + a2 x 4 + a3 x 3 + a4x 2 + a5 x 1 + a6 
dv(x) 



dx 



5 al x 4 + 4 a2 x 3 + 3 a3 x 2 + 2 a4 x + a5 



Las condiciones de frontera de acuerdo a la Figura 4, son: 
x = 0; v(0) = dly x = L/2; v(L/2) = d2y x = L; v(L) = d3y 



x = 0; ^ = OI 

dx 



x = l/2; ^ = 02 



dx 



x = L; *?M = 03 

dx 



Y reemplazando y ordenando en forma matricial se obtiene: 
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(13) 



20 Fernando Olmedo 



Se puede obtener el vector columna de incógnitas [a¡] mediante: 



[K] * [a] = [d] 



[a] = [K]" 1 * [d] 
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El polinomio de interpolación se puede ensamblar multiplicando el vector 
fila O 5 x 4 x 3 x 2 x 1) por el resultado anterior [a] 



/ 



(x 5 X 4 X 3 X 2 X 1) 



24dly 24d3y 4 41 16 42 443 
L 5 L 5 L 4 L 4 L 4 



16 d2y 68 dly 52 d3y 12 41 40 42 8 43 
~ ~ L 4 L 3 L 3 T 3 " 



66 dly 32d2y 34 d3y 13 41 32 42 



L 3 



L 3 
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5 43 



16d2y 23 dly 7 d3y 6 41 8 42 *3 
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/ 



El resultado de esta operación genera un polinomio sumamente grande 
para ser desplegadas en el documento. Este polinomio se reagrupa de tal manera 
que las funciones de interpolación Ni puedan mostrarse en forma explícita y por 
tanto la función de desplazamiento se representara en la siguiente forma: 

v(x) = Ni(x) dly + N 2 (x) efe + N 3 (x) d2y + N 4 (x) 4> 2 + N s (x) d3y + N 6 (x) ^3 (1 5) 

Obteniendo el siguiente resultado para cada función de interpolación N¡: 

Nl(x) = 24 @ 5 -68 0 4 + 66 @' - 23 @ Z + 1 

/X\ 4 /X\ 3 /X\ 2 

N3(x) = 16 (-) - 32 (-J + 16 (-) 

/X\ 5 /X\ 4 /X\ 3 /X\ 2 

NSW = -24 y +52 0 -34 0 +7 (-) 
N2(x) = -4 x Q 4 + 12 x @ 3 - 13 x @' + 6 x @ - x 
N4(x) = -16 x Q 4 + 40 x 0 3 - 32 x 0 2 + 8 x 0 

+ 8x (r) " 5x (r) +x (r) 
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Para L=1 se pueden graficar las funciones de interpolación para verificar 
que tengan el valor de 1 en su nodo asociado y cero en los otros nodos, ver figura 
5 y 6. 



Funciones de Interpolación de alto orden asociadas a bs desplazamientos 




0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figura 5 Funciones de interpolación de alto orden N1 , N3, N5 

Derivas de las funciones de interpolación asociadas a las rotaciones 

0.8 1 ^ . ■ ■ • . 

0.6 • *"7 — — t~ 

0.4 — r ' ' 




' 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

X 

Figura 6 Derivadas de las funciones de interpolación de alto orden N2, N4, N6 



4 OBTENCIÓN DE LA MATRIZ DE RIGIDEZ 



El paso más importante del procedimiento es determinar la matriz de 
rigidez para lo cual se seguirá la formulación dada por el método de Galerkin para 
vigas: 

Donde: K¡j es el término i, j de la matriz de rigidez. 
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Las integrales que deberán desarrollarse son las siguientes puesto que la 
matriz es siempre simétrica. 



Para N 1 

f L/ d 2 Nl (x) d 2 Nl (x) \ 
J o V dx 2 dx 2 ) 

>)dx 
>)dx 
>)dx 
>)dx 

f L/ d 2 N l( x) d 2 N 6 (xA 



* 

dx 2 
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d 2 N!(x) t 
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d 2 N 5 (x) 
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d 2 N 1 (x) + 


d 2 N 5 (x) 
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Para N 2 

r L /d 2 N 2 (x) d 2 N 2 (x) 
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Para N 3 . . . 

f L/ d 2 N 3 (x) d 2 N 3 (xn 

f L/ d 2 N 3 (x) d 2 N 4 (x) \ 

J 0 V dx 2 dx 2 / 

f Lf d 2 N 3 (x) d 2 N 5 (x) \ 



d 2 N 3 (x) d 2 N 5 (x) 



dx 2 



dx 2 



)dx 



Afortunadamente el programa MathCAD posee herramientas de cálculo 
simbólico que permite obtener fácilmente estas integrales, como ejemplo se 
despliega kn y ki2. 
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Obteniendo finalmente la matriz de rigidez y con ella la formulación de 
elemento finito dada por: 
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CARGA DISTRIBUIDA 



Como anteriormente se manifestó, los efectos de las cargas puntuales son 
exactamente evaluados cuando se utiliza la función de desplazamiento cúbica, 
mientras que cuando se utiliza carga distribuida necesariamente se debe 
discretizar la viga, por tal razón el enfoque será con respecto a esta carga. Las 
restricciones y las cargas son aplicadas únicamente en los nodos. La 
aproximación usual es reemplazar la carga distribuida con fuerzas nodales y 
momentos tal que el trabajo mecánico efectuado por las cargas nodales (19) sea 
igual al generado por la carga distribuida (18). El trabajo mecánico debido a la 
carga distribuida w está dado por: 

W distribuida = J 0 L -w * (v(x))dx (1 8) 

Donde v(x) es la función desplazamiento. El trabajo generado por fuerzas 
discretas para tres nodos estaría dado por: 

W discre to = mi d)l + m2 d)2 + m3 d>3 + fl dly + f2 d2y + f3 d3y (1 9) 

Se igualan ambos trabajos según la expresión: 



L 

-w* (NiCx) dly+ N 2 (x) ch + N 3 (x) d2y + N 4 (x) d> 2 + N 5 (x) d3y+ N 6 (x) cj^dx 
= mi d>l + m2 d>2 + m3 d>3 + fl dly + f2 d2y + f3 d3y 



Y se evalúa la integral para cada desplazamiento arbitrario: 

1.- dly = 1, 4)1 = 0, d2y = 0, 4>2 = 0, d3y = 0, 4>3 = 0 

f L 7Lw 
-w * (NiCx) dly + N 2 (x) cfc + N 3 (x) d2y + N 4 (x) c|d 2 + N 5 (x) d3y + N 6 (x) 4> 3 )dx = fl = 

Jo 



30 



2.- dly = 0, 4>1 = 1, d2y = 0, 4>2 = 0, d3y = 0, 4>3 = 0 



f L L 2 w 

-w * (Ni(x) dly + N 2 (x) ^ + N 3 (x) d2y + N 4 (x) 4> 2 + N 5 (x) d3y + N 6 (x) 4> 3 )dx = mi = 

Jo 



60 



3.- dly = 0, 4»1 = 0, d2y = 1, 4>2 = 0, d3y = 0, 4» 3 = 0 



f L 8 L w 

-w * (NiCx) dly + N 2 (x) ^ + N 3 (x) d2y + N 4 (x) 4> 2 + N 5 (x) d3y + N 6 (x) 4> 3 )dx = f 2 = - — 

4.- dly = 0,4>1 = 0,d2y = 0,4>2 = l,d3y = 0,4>3 = 0 

í -w * (N x (x) dly + N 2 (x) ^ + N 3 (x) d2y + N 4 (x) 4> 2 + N 5 (x) d3y + N 6 (x) c|) 3 )dx = m2 = 0 
■¡o 

5.- dly = 0,d)l = 0,d2y= 0,d)2 = 0,d3y= l,d)3 = 0 

f L 7 L w 

-w * (NiCx) dly + N 2 (x) ^ + N 3 (x) d2y + N 4 (x) 4> 2 + N 5 (x) d3y + N 6 (x) 4> 3 )dx = f 3 = - — 
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6.- dly = 0, <j)l = 0, d2y = 0, <j>2 = 0, d3y = 0, <|>3 = 1 
-w * (Ni(x) dly + N 2 (x) + N 3 (x) d2y + N 4 (x) <b 2 + N 5 (x) d3y + N 6 (x) <j> 3 )dx = m3 = - — 

Las fuerzas nodales equivalentes obtenidas están representadas en la 
siguiente figura 7. 



7/w 8lw 7/w 

30 15 30 




(>) (b) 
Figura 7 Fuerzas discretizadas equivalentes 



6 EJEMPLO DE IMPLE MENTACIÓN EN MATHCAD 



Seguidamente se realiza un ejemplo con el fin obtener la solución por el 
método de elementos finitos en base a la nueva matriz obtenida, así como 
comparar la solución tanto con la teoría de vigas como con los resultados 
obtenidos en base del polinomio cúbico. Se realiza por lo tanto el cálculo de una 
viga de acero biempotrada con un solo elemento finito, sometida a una carga 
distribuida w= -100 N/mm, con una longitud L = 1 m, de medidas 1 00 mm x 1 0mm. 

Las propiedades físicas son: E = 200000 MPa, b=10, h=100, I = ^ bh 3 , ver 

figura 8. 




Figura 8 Viga biempotrada sometida a carga distribuida 



Donde las condiciones de contorno son: 

L w 

dly=0 4>1 = 0 m2 = 0 f2y=-8— d3y = 0 tf>3 = 0 
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4 2 2 

L w ^ 4-L-w L w 4-L-w L -w 

384-E I ^5 15 15 15~ , 



Como se puede observar la deflexión en el nodo 2 coincide con el valor 



teórico y max = 



wL 4 
384 El 



para la viga biempotrada. 



Se efectúa la substitución inversa: 




o 
o 

T 4 

L -w 
3&4-E-I 
0 
0 

, 0 




Para determinar las fuerzas y momentos se resta del producto de la 
substitución inversa las fuerzas equivalentes obtenidas en el punto 5. 



{/} = [K){d] - {/„} 



(20) 
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La deformada se la obtiene a través de la función de desplazamiento 
transversal con las condiciones siguientes: 



wL" 
d2y = -384EI' 

x:= 0,0.01.. L 



dly = 0, c()l = 0, d3y = 0, 4)3 = 0 



Nl(x) := 24-| — I -«■ - .66- r - 23- - + 1 



N2(x) := -4-x-| — j + 12.,l z -13.» - ^-x.^-x 



N3(x),l 6 -|- -32. - + 16-- 



N4(.\) := -16-x-l ^) +40-x/^l -32-x/^l + 8-x- 



N5(x) 



-*[ti +52 tij -ní) 3+7 (i 2 
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vFEA(x) = Ni(x) dly + N z (x) + N 3 (x) d2y + N 4 (x) 4> 2 + N s (x) d3y + N 6 (x) 4> 3 (21 ) 



Adicionalmente en el mismo gráfico se dibujara la curva elástica obtenida 
de la ecuación general de vigas 

yED(x)= -J-(-wx 2 )(L-x) 2 (22 ) 

Pudiendo observarse en la figura 9, la total concordancia de ambas 

curvas. 
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Análisis de vigas por el Método de Elementos Finitos 



CURVA ELASTICA, VIGA BIEMP OTRADA CARGA DISTRIBUIDA 



o 

-0.16 
-0.32 
-0.48 
vFEA(x)-0.64 

-0.8 

yED( x) -0.96 
-1.12 
-1.28 
-1.44 
-1.6 



0 



•)()() 



M(x) = E I 



Figura 9 Comparativa curva (21) y (22) 
Igualmente se obtiene el Momento, figura 10 a lo largo del eje mediante: 

d 2 vFEA(x) 



M(x) 



dx 2 



-1x10 



MOMENTO 




100 



20(1 
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400 



500 



0 00 



700 



800 



900 



1x10 



Figura 10 Momento obtenido para la teoría propuesta 

d 3 vFEA(x) 



Y el diagrama de cortante figura 1 1 , con: V(x) = E I 



dx 3 



CORTANTE 



Vcoríant^í) 
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Figura 1 1 Momento obtenido para la teoría propuesta 
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Adicionalmente se efectuará el despliegue de las mismas funciones en la 
figura 12, cuando se utiliza el polinomio cúbico y se dlscretiza la viga mediante 2, 3 
y 4 elementos finitos, utilizando la matriz (12). 



CURVA ELASTICA: 



POLINOMIO CUBICO: 2, 3 y 4 ELEMENTOS 
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vFEA2e(xl) 
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POLINOMIO CUBICO: 2, 3 y 4 ELEMENTOS 
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CORTANTE: 
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POLINOMIO CUBICO: 2. 3 Y 4 ELEMENTOS 
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Figura 12 Curva elástica, Momento y Cortante obtenido con la teoría clásica de elementos 

finitos a partir del polinomio cúbico 

Se observa que cuando se utiliza el polinomio cúbico, no existe 
compatibilidad con la teoría de vigas para el momento y cortante, la compatibilidad 
se logrará únicamente aumentando el número de elementos y así las funciones de 
momento y cortante convergerán con la curva presentada en la figura 1 0 y 1 1 . Con 
4 elementos finitos el valor del momento con el polinomio cúbico es de 
- 7.812 MN mm, con 1 elemento finito el valor del momento usando el polinomio 
propuesto de grado 5 es de -8.333 MN mm. La teoría de vigas predice 



12 12 



-8.333 MNmm 



7 EJEMPLO DE IMPLEMENTACIÓN II 



El segundo ejemplo, figura 13 consistirá en calcular una viga de acero 
biempotrada con un elemento finito, sometida a una carga distribuida w= 
-lOON/mm hasta la mitad de la longitud de la viga, de una longitud L = 1 m, de 
medidas 100 mm x 10 mm. Las propiedades físicas son: E = 200000 MPa, 

b=10mm, h=1 00 mm, I = — bh 3 . Este problema puede realizarse de dos formas: 

utilizando 2 elementos finitos, en ese caso la matriz de rigidez será de 10 x 10 o 
determinar por el método de igualación del trabajo generado por las cargas y 
utilizar un solo elemento finito, se va a optar por este último método. 



B un i í I 

I r 1 

Figura 13 Viga biempotrada sometida a carga distribuida hasta la mitad de la viga 
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Efectuando la Igualación del trabajo mecánico según el punto 5 se 
obtienen las fuerzas equivalentes siguientes según la figura 14: 



r 

1 C 



101/. w 



480 



4lw 
15 



w 

I ! 2| 



33^13w^V 
Z v 960 ^ 



24 



Figura 14 Fuerzas discretizadas equivalentes 
Donde las condiciones de contorno son: 
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Como se puede observar la deflexión en el nodo 2 en este caso es y max = 



wL' 
768 E I 



Se efectúa la substitución inversa: 
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Para determinar las fuerzas y momentos se resta del producto de la 
substitución inversa las fuerzas equivalentes obtenidas en el punto 5. 

{/} = [K]{d] - {/„} 

f 4-L-w N 
15 



15 
8-L-w 
15 
0 

4-L-w 
15 



Finalmente generamos las curvas figura 15 mediante: 

vFEA(x) = N ± (x) dly + N 2 (x) cj^ + N 3 (x) d2y + N 4 (x) $ 2 + N 5 (x) d3y + N 6 (x) 4> 3 
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Igualmente se obtiene el Momento a lo largo del eje mediante: 
M« = EI« 

dx" 1 
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8 CONCLUSIONES 



La utilización del programa MathCAD provee un efectivo y atractivo 
ambiente de aprendizaje del método de los elementos finitos para los estudiantes. 
El método propuesto es de fácil implementación y como se puede apreciar 
sumamente exacto, pudiendo ser aplicado a cualquier problema de ingeniería que 
involucre vigas o ejes. Si bien cuando se utiliza el método clásico el valor de las 
deflexiones no difiere, la incongruencia se da en el cálculo de los momentos y por 
ende de las tensiones por lo que se necesita discretizar. La alternativa a la 
discretización es utilizar elementos de alto orden y la obtención de fuerzas y 
momentos equivalentes de cualquier carga aplicada. 
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